Regressions-Diagnose · 🔍

library(tidyverse)

library(broom)     # tidy model data
library(patchwork) # combine ggplots
library(ggthemes)  # ggplot Stata theme
library(estimatr)  # robust standard-errors

##  tidy lm model results with broom (and ggplot)

tidy_lm <- function(fit) {
  tidy(fit) %>% mutate_if(is.numeric, round, 2)
}

glance_lm <- function(fit) {
  glance(fit) %>% 
    select(r.squared, adj.r.squared, statistic, p.value, df, nobs) %>% 
    mutate_if(is.numeric, round, 2)
}

plot_lm <- function(fit,   bw = 25) {
  dt_pl<-  
    augment(fit) %>% 
    mutate(einkommen_geschaetzt = .fitted,
           residuen = .resid)

  pl1 <- 
    ggplot(dt_pl, aes(x = residuen)) +
    geom_histogram(binwidth = bw, colour = "black", fill = "grey") +
    stat_function(
      fun = function(x) dnorm(x, mean = mean(dt_pl$residuen),
                              sd = sd(dt_pl$residuen))
                        * bw * nrow(dt_pl),
      colour = "blue"
    ) +
   theme_stata()

  pl2 <- 
    ggplot(data = dt_pl, aes(x = einkommen_geschaetzt, y = residuen)) +
    geom_point(alpha = 0.2) +
    geom_smooth(se = FALSE) +
    theme_stata()

  pl1 + pl2 + plot_layout(widths = c(2, 3))
}

Wunsch I · Perfekte Daten

In der Vorlesung 🧑‍🏫 am 19. November 2020 gab es eine Fee 🧚 die Wünsche 🧞 erfüllt.

Den Wunsch 🪄 nach perfekten Daten 📈 erfüllt sie nicht 😥.

Die Fee 🧚 überlegt, nutzt R 👩‍💻 und schafft fast perfekte Daten 📈.

param_const <- 500
param_alter <- 25

n <- 1000
sd_einkommen <-  100

dt_raw <- 
  tibble( 
    id = 1:n,
    alter = seq(20, 60, length.out = n) %>% round(2),
    einkommen_linear = round(param_const + param_alter*alter)
    ) %>% 
  mutate(einkommen = round(einkommen_linear + rnorm(n, sd = sd_einkommen)))

\[ Einkommen_{ID} = 500 + 25 \cdot Alter_{ID} \: + ( \: Fehler_{ID} \: )\]

ggplot(dt_raw, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

Wunsch II · Diagnose für SoWis

Die Studierenden 🙇🙇🙇🙇 hatten einen eigenen Wunsch 🪄 an die Fee 🧚.

Kannst Du uns Regressions-Diagnose 📈 mit einer sozialwissenschaftlichen Geschichte 📚 erklären?

Daten Fee

Wir stellen uns Postbeamte 📯 vor die regelmäßig Lohnerhöhungen 💰 strikt nach dem Senioritätsprinzip 🥸 erhalten.

Wir haben also ein Modell 📈 bei dem das Einkommen 💰 vom Alter 🧑👴 abhängt.

Das wahre Einkommen 💰 (einkommen_linear) ist uns bekannt und unsere beobachteten Daten 🕵️ (einkommen) weichen davon normal-verteilt 🧮 ab

So hat die Fee 🧚 uns die Daten 📊 gegeben.

ggplot(dt_raw, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_jitter(aes(y = einkommen_linear),
              width = 0.1, alpha = 0.5,
              colour = "darkgreen", shape = ".") +
  theme_stata()

DT::datatable(dt_raw)

Modell (fast) perfekt

Wir bekommen schöne Ergebnisse bei der Schätzung des Modells. 😄

fit <- lm(einkommen ~ alter, data = dt_raw)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Annahmen (1-7)

1. korrekte Spezifikation (Linearität)
2. Erwartungswert der Residuen (lokal) = 0
3. Normalverteilte Residuen
4. Varianzhomogenität (gleichmäßige Streuung Residuen)
5. Berücksichtigung einflussreicher Fälle (keine Ausreisser)
6. keine Multikollinearität (geringe Korrelation unabhängiger Variablen)
7. keine Autokorrelation der Residuen (insbes. Zeitreihen)

Korrekte Spezifikation

Dritt-Variablen

Briefträger ✉️ arbeiten oft nur halbtags 🕐 und verdienen 💰 dann auch nur die Hälfte.

dt <-
  dt_raw %>% 
  mutate(einkommen = if_else(id %% 3 == 0, einkommen/2 + rnorm(n, sd = 50), einkommen),
         halbtags = if_else(id %% 3 == 0, "ja", "nein") %>% fct_relevel("ja", after = Inf))

ggplot(dt, aes(x = alter, y = einkommen, colour = halbtags)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Ohne Dritt-Variable 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit, bw=50)

Mit Dritt-Variable 🙂

fit <- lm(einkommen ~ alter + halbtags, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Noch besser mit Dritt-Variable und Interaktions-Effekt 😃

fit <- lm(einkommen ~ alter*halbtags, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Funktionale Form

Bei der Bundespost 📯 wurde eine Bezahlung 💰 nach Humankapital eingeführt.

Es wurde erkannt, dass bei jungen Beamten 🧑 die Produktivität stärker wächst als bei Älteren 🧓. Daher erhalten sie höhere Lohnzuwächse 💰.

dt <-
  dt_raw %>% 
  mutate(einkommen_polynom = -300 + 80*alter - 0.725*alter^2,
         einkommen = einkommen_polynom + rnorm(n, sd = sd_einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
    geom_jitter(aes(y = einkommen_polynom),
                width = 0.1, alpha = 0.5,
                colour = "darkblue", shape = ".") +
  theme_stata()

Ein lineares Modell ist falsch spezifiziert. 😟 Die Regressions-Diagnose zeigt dies. ☝️

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Korrekt ist ein Modell mit einem Polynom (hier einem quadratischem Term). 😃

fit <- lm(einkommen ~ alter + I(alter^2), data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

( Den R Hinweis auf Multi-Kollinearität ignorieren wir einfach. Natürlich gibt es einen Zusammenhang zwischen alter und dem quadrierten alter^2. )

Varianzhomogenität

Bei der Post 📯 haben hohe Beamte und Minster ein höheres Einkommen 💰 als einfache Beamte. Aber nicht alle werden hohe Beamte 🧑‍💼 oder Post-Minister 🏤.

dt <- dt_raw %>% mutate(einkommen = einkommen_linear + id/500 * rnorm(n, sd = sd_einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
    geom_jitter(aes(y = einkommen_linear),
                width = 0.1, alpha = 0.5,
                colour = "darkgreen", shape = ".") +
  theme_stata()

Im Standard-Modell erhalten wir verzerrte Standard-Fehler. 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Besser ist ein Modell mit robusten 💪 Standard-Fehlern. 😄

fit <- lm_robust(einkommen ~ alter, data = dt)
tidy_lm(fit)

glance(fit) %>% mutate_if(is.numeric, round, 2) 

Multikollinearität

Ältere 👴👵 Postbeamte 📯 sind kleiner als Jüngere 👧👦.

dt <- dt_raw %>% mutate(groesse_cm = round(seq(180, 160, length.out = n) + rnorm(n, sd = 2.5), 3))

ggplot(dt, aes(x = alter, y = groesse_cm)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

Es gibt dann auch einen Zusammenhang 📉 zwischen groesse_cm 👵 und einkommen 💰, natürlich keinen kausalen ⚙️.

ggplot(dt, aes(x = groesse_cm, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

fit <- lm(einkommen ~ groesse_cm, data = dt)
tidy_lm(fit)

glance_lm(fit)

Nehmen wir alter 👵 und groesse_cm 📐 in einem Modell 📉 auf haben wir Multikollinearität 🔗 😟

fit <- lm(einkommen ~ alter + groesse_cm, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Multikollinearität im Modell 📉 erkennen 🕵️ wir auch an der sehr starken Korrellation der unabhängigen Variablen.

dt %>% 
  select(alter, einkommen) %>% 
  cor() %>% 
  as_tibble() %>% 
  mutate_all(round, 2)

Ausreisser

Bei der Datenübertragung 📡 gab es einen Fehler ⚠️ und die letzte Ziffer 🔢 wurde bei einigen Gehaltsangaben 💰 entfernt.

id_random <- runif(80, min = 1, max = n) %>% round()

dt <-
  dt_raw %>% 
  mutate(einkommen = if_else(id %in% id_random, round(einkommen/10), einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Wir erhalten eine verzerrte Schätzung. 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit, bw=50)

Die falschen Beobachtungen müssen wir entfernen (oder korrigieren). 😃

dt_tmp <- dt %>% filter(einkommen > 500)

fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Autokorrelation

Die Postbeamten 📯 wurden 2010 und 2020 📅 nach ihrem Einkommen 💰 befragt. ☎️

dt_tmp <-
  dt_raw %>% 
  mutate(einkommen = round(einkommen - param_alter*10 + rnorm(n, sd = 10)),
         alter = alter - 10,
         befragung = 2010) %>% 
  filter(alter >= 20)

dt <-
  dt_raw %>% 
  mutate(befragung = 2020) %>% 
  rbind(dt_tmp) %>% 
  arrange(id, befragung) %>% 
  relocate(befragung, .after = id)

DT::datatable(dt %>% filter(id > 500))  # Daten ab ID 500 um Panel-Struktur zu verdeutlichen 

Postbeamte 📯 die 2010 📅 älter als 50 waren sind inzwischen in Pension 👵 und haben an der 2020 Befragung ☎️ nicht mehr teilgenommen.

dt_tmp <- dt %>% mutate(befragung = factor(befragung))

ggplot(dt_tmp, aes(x = alter, y = einkommen, colour = befragung)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Im Modell (n = 1750) sind Beobachtungen (und Residuen) nicht mehr unabhängig. 😟 Wir haben ja viele zweimal befragt. ☎️

fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)

glance_lm(fit)

plot_lm(fit)

Aufgabe

Diskutieren Sie 🙇🙇🙇 die Regressions-Annahmen 📈 und zeigen Sie ob und wie Verletzungen 🏥 der Annahmen erkannt 🕵️ werden können.

---

Quellen

• https://commons.wikimedia.org/wiki/File:SophieAndersonTakethefairfaceofWoman.jpg
• https://commons.wikimedia.org/wiki/File:Bundesarchiv_Bild_183-63142-0001,_Cottbus,_Brieftr%C3%A4gerinnen.jpg
• https://commons.wikimedia.org/wiki/File:Bundesarchiv_Bild_183-1990-0924-020,_Gera,_Protest_von_Post-Gewerkschaftern.jpg
• https://commons.wikimedia.org/wiki/File:Bundesarchiv_Bild_183-U0205-0010,_Berlin,_Tag_des_Post-_und_Fernmeldewesens,_Vorbereitung.jpg
• https://commons.wikimedia.org/wiki/File:Bundesarchiv_Bild_183-59111-0002,_Berlin,_Postamt_N_58,_Brieftr%C3%A4ger.jpg
• https://commons.wikimedia.org/wiki/File:Telefonzentrale_bei_der_2.Infanteriedivision_(BildID_15705125).jpg
• https://commons.wikimedia.org/wiki/File:Radioquiz_%22Allein_gegen_alle%22_aus_dem_Ratssaal_(Kiel_45.567).jpg

---